%\documentclass{llncs}
%\documentclass{article}
\documentclass[12pt,a4,version]{article}
\usepackage[margin=0.85in]{geometry}

\usepackage{amsmath,amssymb,latexsym,graphicx,rotating,algorithmic,algorithm}
\usepackage{url}
%\usepackage{tikz}
\usepackage{amssymb}
\usepackage{graphics,color}
%\usetikzlibrary{snakes}
\def\BO{\mathcal{O}}
\def\NP{\mathcal{NP}}
\def\LL{\mathcal{L}}
\newcommand{\cC}{\mathcal C}
%\newcommand{\comment}[1]{\textcolor{blue}{\textbf{#1}}}
\newcommand{\shlomo}[1] {\textcolor{red}{shlomo: #1}}
\newcommand{\yossi}[1] {\textcolor{blue}{yossi: #1}}
\newcommand{\igor}[1] {\textcolor{green}{Igor: #1}}

\newtheorem{theorem}{Theorem}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{defn}[theorem]{Definition}
\newtheorem{eg}[theorem]{Example}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{rem}[theorem]{Remark}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{construction}[theorem]{Construction}

\newenvironment{corollary*}
{\vspace{0.5 em} \noindent {\bf Corollary} \em}{\vspace{0.5 em}}

\newenvironment{definition}{\defn \em }

\newenvironment{example}{\eg \em }

\newenvironment{remark}{\rem \em }

\newenvironment{proof}
{\vspace{0.5 em} \noindent {\bf Proof \,}}{\QED {\vspace{0.5 em}}}

\def\filledbox{\vrule height 1.8ex width .8ex depth 0ex }
\newcommand{\qed}{\filledbox}

\def\QED{\ifmmode\squareforqed\else{\unskip\nobreak\hfil
\penalty50\hskip1em\null\nobreak\hfil\qed
\parfillskip=0pt\finalhyphendemerits=0\endgraf}\fi}

\newcommand{\ignore}[1] { }


\begin{document}
\date{}



\title{Mixed Integer Programming Accelerator\\
{\normalsize M.Sc. Research Proposal} }
\vspace{5mm}

\author{Igor Zarivach \\ \\ \\
Advisors: \\
\\
Prof. Shmuel Zaks \\
Technion  \\
 Israel Institute of Technology \\ \\ \\
Prof. Shlomo Moran \\
Technion  \\
 Israel Institute of Technology \\ \\ \\
\\
Dr. Yossi Shiloach \\
\\
\\ \\ \\
\\ \\
\today
}
\maketitle

\vspace{5mm}
\thispagestyle{empty}
\clearpage

\tableofcontents

\thispagestyle{empty}
\clearpage
\setcounter{page}{1}
%\vspace{5mm}
%\vspace{-1 em}
\section{Introduction} \label{sec:intro}
Optimization is the task of finding an optimal solution to a problem that has many alternative solutions. The feasibility is defined with respect to a large set of constraints that any solution must satisfy.
The optimality is defined with respect to a certain measure that can be assigned to each feasible solution. \\
Mathematical optimization deals with the mathematical representation of an optimization problem. Transforming the field problem into a mathematical representation is called ‘modeling’ and is an art by itself. The model has to be logically equivalent to the original problem and express the user relative priorities as closely as possible in the objective function. It also has to support efficient search.
Different models to the same problem may be both correct and logically equivalent but may converge into an optimal solution in totally different times.

There are a handful of optimizers or solvers, computer programs that solve optimization problems, such us IBM CPLEX, Gurobi and the others. The user models the problem using the modeling language and supplies the data. The solver uses mathematical algorithms to find the optimal solution.

Optimization problem with integrality constraints is called a Mixed Integer Programming (MIP). Solving MIP for optimality is NP-hard ~\cite{garey2002computers}.  Recent advances in the implementation of algorithms and hardware improvements have extended the size of MIP programs that can be solved efficiently using today's computers.

However, there are still hard problems that can't be solved efficiently. MipLib, an open library of MIP problems, contains a standard test set used to compare the performance of MIP optimizers.
MIPLib classifies instances into three categories: easy for those that can be solved
within an hour on a contemporary PC with a state of the art solver, hard for those
that are solvable but take a longer time, and finally open instances for which the optimal solution is not known. As for 2010, there were still 42 hard, and 134 open problems. \cite{MipLib}

When the problem is computationally hard, we are looking for the best solution possible. This research proposes a new approach or heuristic for finding "good" feasible solutions solving general MIP models.


\section {Linear programming}
The most frequently used method of mathematical optimization is called linear programming (LP). Only linear constraints are allowed as well as linear objective function. It turns out that
many optimization problems can be modelled by linear constraints.

Below is a mathematical formulation of linear minimization problem of $n$ variables and $m$ constraints using matrix notation:
\begin{align*}
\text{min} \quad
c^T x& \\
\text{s.t.} \quad
Ax &\le b \\
x &\geq 0
\end{align*}
where $x=(x_1,..,x_n)$ are variables, $c=(c_1,...,c_n)$ are the coefficients of the objective function, A is $mn$ matrix and $b=(b_1,..,b_m)$ are constants.

The linearity of the mathematical representation allows very efficient algorithms for searching the
optimal solution. The best known is Simplex designed by Danzig  ~\cite{dantzig1951maximization}. Simplex method proceeds by
walking from one vertex to another of the convex polytop defined by linear inequalities. At each step, it walks to a vertex that improves the objective function and finally arrives to the global optimum or decides that the problem is not feasible.\\
Another linear programming algorithm is the interior-point method introduced by Karmarkar ~\cite{karmarkar1984new}. Interior point method is polynomial and its implementations are computationally competitive with simplex in practice.

\section{Mixed integer linear programming}
If variables are restricted to take integer values, the problem is called \emph{integer linear programming}. When some variables are integer and some are allowed to take real values, then the problem is called a \emph{mixed integer linear programming} (MILP or MIP).\\
\emph{Linear programming relaxation} of a MIP problem is achieved by removing integrality constraints. The resulting problem is LP.

MIP is one of the most largely used optimization method, the modeling is used in production scheduling, transportation, logistics, bioinformatics and other areas ~\cite{Williams}.
Whereas the simplex method is effective for solving linear programs, there is no single technique for solving integer programs efficiently.

Current state of the art MIP solvers can't solve many real world problems to optimality. However, feasible solutions with good objective value can still provide a significant practical value, for example:
\begin{enumerate}
\item Practical problems must be solved within time constraints which don't allow to solve for optimality.
\item Feasible solution objective value is satisfactory.
\item Feasible solution is close enough to the optimal.
\end{enumerate}
It is therefore often important to find the best possible solution as fast as possible. There are different heuristics deployed by most solvers that find feasible solutions with "good" enough objective value. Benchmarking helps to compare various strategies.

\section{MIP Solving Methods}
 The MIP optimization problem is NP-hard ~\cite{garey2002computers}, so there is no known polynomial algorithm that finds the optimal integer solution. However, there are methods that can find near optimal integer solution effectively in practice.

\subsection{Cutting planes}
This method works by adding additional constraints to the model. The MIP problem is relaxed and the relaxed LP problem is solved. If the solution is fractional, we add a new constraint which removes the fractional optimal solution, but doesn't change the feasible set of integral solutions.
Then, the process is continued until an integer solution is finally found.
Cutting plane methods have proved not very successfull on large problems.

\subsection{Branch and bound}
Branch and bound is a divide and conquer technique for finding the optimal solution on a general MIP problem ~\cite{Wolsey}. Lets consider a minimization optimization problem.
Branch and bound runs LP solver on the original MIP problem without integrality constraints (LP relaxation).
If the solution is integral and feasible for the MIP problem, this is the optimal one. Otherwise, the idea is to divide the problem to smaller MIP problems (branching) and optimize them recursively.\\
The recursion stops on sub-MIP problems that should not be divided further. Main ideas of the branch and bound are summarized below.

\paragraph{The tree}
The solver enumerates solutions using the enumeration tree. Nodes of the tree represent MIP problems. The tree grows from the root, the original problem with all feasible solutions set. Every time a problem is divided to new MIP subproblems, new child nodes are added to the tree. Child nodes represent the MIP problem of the parent node with an additional constraint.

\paragraph{The incumbent solution and pruning}
The incumbent solution is the best feasible solution found so far in the process of branch and bound. Branch and bound finds feasible solutions either from solving LP relaxation or by heuristics.
Having the incumbent solution with good objective value will improve the search performance significantly by discarding subproblems which have inferior solutions to the incumbent. Pruning is one of the most important aspects of branch and bound since it is precisely what prevents the search tree from growing exponentially.

\paragraph{The search algorithm}
The solver traverses the tree starting from the root. The traversal decides whether the node should be branched or pruned. While at some tree node, the solver uses LP relaxation to find a lower bound on the optimal solution in the current sub MIP problem. There are few cases:
\begin{itemize}
\item The problem is unbounded. Then the original problem is unbounded and the solver stops the traversal.
\item The problem's optimal solution is not better than the best solution found so far, the incumbent.
Then there is no need to continue the tree traversal from the current node. The current node is pruned.
\item The problem's lower bound improve the best solution found so far. Then the traversal continues from the current node and branching will be later performed.
\item The problem is infeasible. The node is pruned.
\end{itemize}

\paragraph{Branching}
Given a fractional solution to LP relaxation, the algorithm chooses some variable $x_i$ that is restricted to be integer, but whose value in the LP relaxation is fractional $f$. Then the problem is divided to two different subproblems:
\begin{enumerate}
\item By adding a constraint $x_i\le \lfloor{f}\rfloor$
\item By adding a constraint $x_i\ge \lceil{f}\rceil$
\end{enumerate}

\paragraph{Lower and upper bounds gap}
During the search, the solver maintains both lower and upper bounds on the optimal solution.  The lower bound is obtained by taking the minimum of the lower bounds of the current leaf nodes. The upper bound is the lowest objective value solution found so far, the value of the incumbent.
The difference between the current upper and lower bounds is known as the gap. When the gap is zero, the solver has found the optimal solution.

\section{Solvers}
There are many MIP solvers, both commercial such as IBM ILOG CPLEX, Gurobi and open source. We have chosen IBM ILOG CPLEX for its state of the art solver and free usage for academic research purposes.

\subsection {IBM ILOG CPLEX}
IBM ILOG CPLEX is the commercial solver by IBM. It is part of the IBM ILOG Optimization Suite of tools, which contains also the modeling language (OPL).

\paragraph{OPL}
OPL is a declarative language for defining mathematical models. With OPL the user declares data elements, decision variables, objectives and constraints. OPL supports pre and postprocessing, as well as flow control using IBM ILOG Script.

\paragraph{CPLEX}
CPLEX includes many solvers, such as the LP simplex and interior point, MIP optimizer (branch and bound) and concurrent optimizer which runs several algorithms parallelly.

\paragraph{OPL APIs}
IBM ILOG OPL Interfaces provide the developer with programmatic access to OPL model and CPLEX solver from external applications. Using OPL API's it is possible to control the CPLEX solver from external Java, C++ or .NET application.


\section{Research goal}
Given a general MIP problem, we want to find a feasible solution with good objective value in a given amount of time. The feasible solution objective value should be superior to the solution found by state of the art IBM CPLEX solver in the same time period.
The quality of the solution will be compared using the benchmarking on several models against CPLEX solver.

\section{Related Work}
Most commercial solvers use heuristics to find "good" solutions in a short time. The heuristic process is integrated into the branch and bound tree search.

\paragraph{Rounding and diving} During the branch and bound tree search, the solver finds a fractional solution to LP relaxation problem and tries to convert it into the integer feasible solution. This is obtained through rounding an LP solution either component by component or by \emph{diving}, rounding one or more variables and solving the LP relaxation again, iteratively.

\paragraph{Local search}
Once one or more feasible solutions are available, local search algorithm improves its objective value.
The algorithm defines a neighborhood of the solution, and replaces the solution by its neighbour with improved objective value. The process stops when the algorithm gets a locally optimal solution.

\paragraph{Solution polishing}
Solution polishing, described by Rothberg ~\cite{rothberg2007evolutionary}, explores the neighborhood of given feasible solution by fixing some of the binary variables, and solving a sub-MIP problem.
Polishing uses a suite of evolutionary algorithms, combined with large-neighborhood search to explore
possible mutations and combinations of a given solution population. The approach was proved to find significantly better solutions than the best existing heuristics for several very difficult MIP models ~\cite{rothberg2007evolutionary}. \\
CPLEX implements solution polishing in a second phase after a first phase of conventional branch and cut, due to high cost of calculations. Only after some feasible solution has been found, solution polishing can improve the best known solution.

\section{Research plan}
We intend to implement the MIP accelerator software that finds near optimal MIP solutions \emph{faster} than off the shelf MIP solvers. The software will search for a good feasible solution by solving many smaller sub MIP problems with a standard MIP solver. \\
Given a MIP problem and some initial feasible solution, we create a new sub MIP problem by fixing a set of integer variables to their values in the feasible solution. The objective value and constraints are the same as in the original problem. This sub-MIP problem has smaller number of integer variables, thus the size of the branch and bound tree is reduced and its search can be done efficiently by MIP solver.  \\
Running the MIP solver on the sub MIP problem produces a feasible solution to the original MIP problem. Because the tree size is smaller, the MIP solver should find a feasible solution with better objective value than the initial one or decide that none exists. Having a better feasible solution will improve the branch and bound search time in the next iteration by effective pruning of nodes.
Repeating the process with various fixing strategies, we can improve the objective value of feasible solution fast.
This fix and run process is different from exisiting heuristics that are applied by the solver during the branch and bound tree search. There are many sequential branch and bound searches initiated by the MIP accelerator for each of the sub MIP problems.
The accelerator will run in iterations, improving the feasible solution.
Each iteration runs the following algorithm:
\begin{algorithm}[H]
\caption{FixAndImprove ($M,f$) - improve objective value of the given feasible solution}
\label{alg1}
\begin{algorithmic}[1]
\REQUIRE MIP problem $M$ = min $c^Tx$ : $Ax \le b$, $x_i \ge 0$, $x_i \in \mathbb{Z}$
\REQUIRE A feasible solution $f=(f_1,...f_n)$ to MIP $M$
\ENSURE {Return a feasible solution $f'$ with better objective value than $f$, $c^Tf'$ $<$ $c^Tf$}
\STATE $I\leftarrow$  Choose a subset of integer variables to be fixed in $f$, $I\subseteq\{1,2,...,n\}$
\STATE Generate a new problem $M'$: add to $M$ equalities $x_i=f_i$, $\forall i\in I$
\STATE Call the solver to solve $M'$ for predefined time period
\STATE $f'\leftarrow$ solution for $M'$
\RETURN $f'$
\end{algorithmic}
\end{algorithm}

\subsection{Fixing strategy}
The amount and identity of  variables to be fixed affects the possible extent of the improvement and the computation time of each iteration. When more variables are fixed, solution space gets smaller. This allows fast implicit enumeration by branch and bound search. However, the quality of optimal solution might be inferior. On the other hand, small amount of fixing will require long enumeration process by branch and bound.\\

Which variables are suitable for fixing? The algorithm knows which variables are integral from the model.
However, not all variables are beneficial to fix. For example, variables that are computed from other variables should not be fixed.\\
We propose two strategies for the fixing:

\paragraph{Fix randomly}
Choose variables randomly from the set of all integer variables. This is the straightforward implementation of the fixing. A fraction of variables to be fixed will be provided as the parameter to \emph{FixAndImprove} [Algorithm  \ref{alg1}].

\paragraph{Fix variables based on domain knowledge}
We can use the information provided in the model (\emph{domain knowledge}) to decide which variables are candidates for fixing.
Actually, the user will define a group of variables that are related in the problem domain. \\
For instance, in scheduling problem of jobs to machines, a binary variable $x_{ij}$ represents a decision to assign job $j$ to machine $i$. All variables $x_{kj},\forall j$ for a given $k$ are related since they represent assignments to $k$-th machine. So, the user can define a group of variables $M_k = \{x_{kj}|\quad \forall j\}$ using the modeling language.
Given definitions of groups of variables, the accelerator can have better fixing decisions using domain knowledge insights. One option is to fix every variable in the group. \\
For the scheduling jobs problem, this reduces the problem size by one already assigned machine.\\

How many variables do we fix in each iteration? We can decide based on the quality of the feasible solution. During first iterations, a small fixing fraction can produce good feasible solutions because most of the solution space is unexplored. If after several iterations there is no improvement in the feasible solution objective value, we should reduce the fixing fraction and increase the solution space.

\subsection {Attributing}
Problem attributing is a preprocessing phase in which the user assigns certain attributes to variables. These attributes bear problem domain knowledge and convey it to the accelerator.
Using these attributes, the user will define groups of variables to be fixed. There are two options for implementing attributing:
\begin{enumerate}
\item Using standard modeling language, such as OPL. The accelerator will use OPL APIs to query the model for attributes.
\item Using proprietary data format. The accelerator will read and parse the data.
\end{enumerate}

\subsection {Solution pooling}
The accelerator needs some feasible integer solution for the fix and search process. Usually, the best feasible solution from the previous iteration is chosen as the input for the next iteration.
However, this approach might restrict the neighbourhood of the search and arrive at local but not global optimum.
Another option is to maintain a pool of feasible solutions during the run of accelerator. Then we could run the fix and search process from any of them and choose the best solution found.


\subsection{Software development}
The software will be developed in Java programming language using agile development process.
The following modules will be developed:
\begin{enumerate}
\item Fixing strategy. Will be responsible for the decision of which variables are to be fixed.
\item CPLEX interfaces. Will provide an access to CPLEX API's.
\item Model. Will be responsible for interaction with a mathematical model.
\item Solver controller. Will be responsible for executing the main search algorithm, by running CPLEX solver iteratively on modified models.
\item View. Will be responsible for interaction with the user and displaying status and results.
\end{enumerate}

\subsection{Applications to Distance based Phylogenetic Reconstruction}
Phylogenetic analysis reveals the relationship among different organisms and their evolutionary history.   Phylogenetic tree represents the evolutionary events that created existing organisms: the leaves of the tree are living species, while internal nodes represent hypothetical common ancestors.
Distance based reconstruction methods attempt to reconstruct the evolutionary tree from the pairwise {\em evolutionary distances} between the species, which are usually estimated by the differences between aligned DNA sequences of the species, (for example, see Chapter 11 in  ~\cite{felsenstein2004inferring}).\\
Although there are specialized polynomial algorithms for distance based tree reconstruction, such as neighbor joining (NJ)~\cite{saitou1987neighbor}, solving the problem using MIP has several advantages:
\begin{enumerate}
\item NJ focuses on the solution (the algorithm). MIP focuses on the problem (the model).
\item NJ is harder to implement and modify, comparing to the compact and easily modifiable MIP models.
\item MIP enables to give different treatment to different distances: An inherent difficulty in distance based reconstructions is that the errors in distance estimations grow exponentially with the distance. Hence most reconstruction algorithms simply ignore distances that are too large. This significantly complicates the algorithm, since the common reconstruction algorithms like Neighbor Joining assume that all pairwise distance estimations are given. \\
On the other hand, MIP enables to easily extend this binary approach to a continuous approach in which the penalties for deviations from the estimated distances decreases as the distance estimation grows. In addition, unlike the common neighbor joining algorithms, MIP models are easily adjusted to handle missing distances.

\end{enumerate}
We will use the phylogenetic MIP model as a benchmarking instance for performance comparison of the accelerator.

\subsection{Benchmarking}
We will compare the results of running the accelerator for a given time versus CPLEX with Polish heuristics on several MIP models. The models will be taken from the MipLib instance library.\\
We also plan to run the accelerator on the MIP model of evolutionary distance based tree. The distance data will be created using a simulation of tree evolution with a probabilistic model. Again, the objective value of the best tree found by the accelerator will be compared vs the CPLEX solver.\\
Another test is to compare the performance of the \emph{fix randomly} and \emph{fix based on domain knowledge} strategies.

\nocite{*}
\bibliographystyle{alpha}
\bibliography{IgorThesisProposal}


%v0.2
\end{document}
